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SECTION I 
INTRODUCTION 


1.0 BACKGROUND 
1.1 Orientation 

Because of their inherent flexibility of control, reliability, and 
efficiency, uses of electromechanical actuators and electric drives are 
finding Increased Interests in aircraft applications. Three areas of 
technology advancement over the past decade are largely responsible for 
placing electromechanical energy converters into a favorable position 
when compared to hydraulic motors or actuators: 


(a). 

Development of rare earth permanent magnet (FM) motors with 
Inherent high efficiency and high power to weight ratios. 

(b). 

Advancements in power level solid-state devices with high 
switching speeds. 

(c). 

Emergence of microprocessors allowing control capability 
and sophistication far surpassing that of analog systems 
while reducing the volime occupied. 

The conventional or commutator DC machine performance characteristics in 
the areas of speed control and position control are highly desirable. 

The brushless DC machine has wide range torque-speed characteristics like 
unto the commutator DC machine without the commitator-brush maintenance 
problems. In addition, the brushless DC machine with a permanent magnet 
rotor has certain superior features to the conventional DC machine: 


Field excitation is eliminated, which removes the complexity 
of supplying power to a rotating meiid>er. Also, machine 
efficiency is increased due to absence of field excitation 
losses. 

(b). 

Higher speed design is possible for FM rotors than is feasible 
with WDimd rotors permitting Increased gear ratios, which leads 
to substantial reduction in electric machine power to weight 
ratios. 

(c). 

Thermal transfer characteristics are improved since the bulk 
of the losses (ohmic and core losses) are generated within the 
stationary member, allowing efficient implementation of fluid 
cooling. 
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1.2 Identification of Problem 

Brtishless DC motor performance is reported in the literature, but it 
concentrates on the nature of Instantaneous voltage and current along with 
average values of developed torque Cl-93. These works typically describe 
systems that do not operate continuously at near zero speed, and thus, 
only average value of torque is of concern. However, the brushless DC 
machine Inherently has oscillatory components of Instantaneous torque 
at all frequencies that are Integer multiples of six times the electrical 
radian of the stator Impressed voltage. Since stator frequency Is directly 
related to rotor position, these oscillatory components are In the fre¬ 
quency range of mechanical system response at low speed values. 

Published works relating to performance of brushless DC machines 
have not analyzed cases of position or sustained near zero speed operation; 
thus, pulsating torque components have been considered of no significant 
consequence. However, Williamson et al CIO3 have mentioned that torque 
capability deteriorates at low speeds, and Demerdash and Nehl have shown 
some Instantaneous developed torque and power wave forms Cll3 without 
comment on the oscillatory components. 

Widespread usages of brtishless DC motor drive systems In low-speed 
and position-control actuator applications are contingent upon development 
of control philosophies and hardware realizations of power conditioning 
arrangements that allow bidirectional power flow while resulting In 
instantaneous motor developed torques that are free of harmonics in the 
range of response for coupled mechanical loads. 

1.3 Consideration of Alternative Technologies 

There are two basic power electronic circuits that are used as power 
conditioning links to couple brushless dc motor drives to a high frequency, 
three-phase AC source: 

(a) DC link Inverter 

(b) Cycloconverter link 

Power circuits of these two basic approaches are Illustrated In functional 
block form by Figure 1. Variations of each arrangement are made .according 
to the control philosophy adopted to satisfy required motor performance. 
Further, for low voltage (500 V or less), low current (200 A or less) 
applications, power circuits can be synthesized with transistors for 
controlled switching elements, while for high power level applications 
silicon controlled rectifiers (SCRs) must be used for controlled switching 
elements. For this study, SCR switching elements are assumed. 

1.3.1 Operation for Speeds Not Hear Zero . For the case of a DC link 
Inverter driving a brushless DC motor that does not operate at near zero 
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Figure 1. Alternative Power Circuits 

(a) DC Link Inverter 

(b) Cycloconverter Link 
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speed, the rectifier may be a diode bridge while the chopper Is used to 
vary magnitude of DC voltage applied to the inverter Input terminals; or, 
the chopper may be eliminated and a phase-controlled converter used as 
the rectifier to vary Inverter input voltage. 

For the case of a cycloconverter, the necessity of a DC link does 
not exist, and thus, there Is only one power conditioning module Inter¬ 
connecting the three-phase AC source with the brushless DC motor. 

Generally, the cycloconverter may be realized as either a midpoint or a 
full-bridge arrangement and either phase-control or synchronous envelope 
operation Implemented. 

1.3.2 Operation for Speeds Near Zero . For near zero speed operation of 
a DC link Inverter for high power level systems, the counter emf of each 
motor phase Is small enough so that natural commutation of the Inverter 
SCRs Is prohibited. In such a case, the chopper Is used to reduce 
Inverter terminal voltage to zero sufficiently long to accomplish commu¬ 
tation of the Inverter SCRs C3]. Since the chopper must be controlled 
so as to enhance Inverter commutation, the rectifier must be a phase- 
controlled converter to permit necessary Inverter Input voltage magnitude 
control. 

For near zero speed operation of a phase-controlled cycloconverter 
tislng circulating current free mode, discontinuous load current tends to 
occur, leading to Increase In load current harmonics. Continuous 
current can be restored by changing to circulating current mode of control 
at the expense of Increased losses due to circulating reactive current Cl23. 
Control Is generally simplified at fractional hertz frequency operation by 
use of synchronous envelope control with circulating current free mode. 

1.3.3 Operation for Regenerative Power Flow . For the case of a DC link 
Inverter, regenerative power flow requires further power circuit modifica¬ 
tion. Gating signals to Inverter SCRs are suppressed and the associated 
shunting diodes function to form a full-bridge diode rectifier. A 
switching circuit Is Introduced in the DC link to reverse polarity of the 
voltage appearing at the DC terminals of the rectifier module. The 
rectifier must be a phase-controlled converter operated In the synchronous 
Inversion mode. Reversal of motor developed torque requires three 
coordinated control actions: phase forward and suppression of Inverter 
SCRs; polarity reversal switch activation; and, phase forward of phase- 
controlled converter SCRs. 

Regenerative power flow using the cycloconverter link Is Introduced 
by simply delaying SCR firing angles beyond 90°. Rate of change of motor 
developed torqvie from a positive to a negative value Is controlled by the 
rate at which the SCR firing angles are changed. Thus, the nature of 
transition from motoring to regeneration Is determined by a single control 
action leading to smooth change with minimum time delay. 
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2.0 OBJECTIVE 


Of the presently available power conditioning technologies for 
linking the brushless DC motor to a high frequency polyphase AC source in 
a near zero speed or actuator application, the cycloconverter using syn¬ 
chronous envelope control with circulating free current mode of operation 
appears most promising In that it requires no external commutation cir¬ 
cuitry and it can be smoothly and quickly changed from motoring to 
regeneration by a single control action. This arrangement using a mid¬ 
point, three-pulse cycloconverter is chosen for study with pulse width 
modulation (FWM) control to eliminate those harmonics from the output 
waveform that lead to harmonics In the developed torque In the range to 
which coupled mechanical loads can respond. 

3.0 SCOPE 

The following tasks were carried out to accomplish the objective: 

(a) . Determine an appropriate mathematical model of the system for 

near zero speed operation. 

(b) . Determine the current harmonics to be eliminated In order to 

remove undesireable pulsations In motor developed torque. 

(c) . Develop a modulation function to use In control of phase 

voltage for eliminating the undeslred current harmonics. 

(d) . Calculate motor performance over a wide speed range with and 

without the FWM phase voltage. 

(e) . Compare ohmic losses with and without FWH to assess the effect 

of harmonic elimination on efficiency. 

















SECTION II 


PROCEDURE 


1.0 SYSTEM DESCRIPTION 

Typical values for motor parameters were obtained by ratio of known 
valxies reported In the literature Cll] where the results are given In 
Table 1. The motor Is wye connected with a maximum design speed of 
45000 rpm. 


TABLE 1. MOTOR PARAMETERS 
Parameter Symbol Value 


No. poles 
EMF constant 
Stator resistance 
Stator leakage Inductance 


P ^ 

K 0.0225 V-s/rad 

R. 0.4 0 

La 25 uH 


Characteristics of the balanced three-phase source are listed In 
Table 2. 


TABLE 2. SOURCE CHARACTERISTICS 


Quantity 


Value 


Phases 3 

Line voltage 138 V 

Frequency 7950 Hz 


The power circuit schematic of a midpoint, three-pulse cycloconverter 
Unking the high frequency, three-phase AC source to the brushless DC 
motor Is depicted by Figure 2. Use of the switch In the midpoint or 
neutral line will be discussed later. 

2.0 ASSUMPTIONS 

When analyzing rare earth permanent magnet machines with non-magnetlc 
retaining rings. It has been found that position dependence and Interphase 
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coupling of machine Inductances can be neglected [11]. Such an approxima¬ 
tion leads to decoupled equations that can be used in networks formed by 
addition of power conditioning circuitry with minimum difficulty. Motor 
counter emfs are taken to be sinusoidal. Further, since the study of 
this report Is concerned with low speed operation, the commutation overlap 
Is very small compared with a period of the motor current wave form and 
can be neglected. 

The SCRs are modelled by a 0.02 ohm forward resistance (2 V at 100 A) 
and a blocking resistance described by 1000 1^ ohms. 

3.0 MATHEMATICAL MODEL 

Referring to Figure 2, using the above assumptions, and applying 
Klrchhoff's voltage law results In the following set of simultaneous 
differential equations to describe the motor electrical performance: 


dt 


r <-‘ 3 “ * ’'cn - *3) 


(3) 


Ri, R 2 , and R 3 are the sum of the motor phase resistance and the resistance 
of the partlcxilar SCR between the source and phase of the motor at the 
particular Instant of time under analysis. The form of phase voltages 
(^an» ^bn» ^cn^ discussed later. 

Neutral connection current Is given by Klrchhoff's current law as 


inN - 

As long as motor speed (u^) remains non zero, the electromechanical 
developed torqiie Is given by 


Tj - (ejij + e^^ + ejij)/^^ (5) 

Equations (1) - (5) completely describe the Instantaneous performance 
of the brushless DC motor. 
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4.0 CURREirr HARMONICS TO BE ELIMINATED 

For both the case of a DC link Inverter drive and the case of a 
synchronous envelope cycloconverter drive, the motor phase currents 
approach the pulse width controlled square wave of Figure 3. A Fourier 
series representation of the wave form shows that all odd, nontrlplen 
harmonics exist In the phase current. 

1, ■ — y 1. cos (HI.) sin nut, n - 1,5,7,11,13, •• • (6) 

1 ^ n “ 6 



Figure 3. Approximate Form of Motor Phase Current 

Similarly, the remaining balanced phase currents are given by 

1^ - M y i cos (Hp sin n(ut - -i^) (7) 

I 3 - ^ y ^ cos (Hp sin n(ut + ( 8 ) 

IT n “ o 3 

n - 1,5,7,11,13,••• 

Since the motor counter emfs sinusoids of fundamental 

frequency and form a balanced three-phase set, each can be expressed as 
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Ejj sin («t + 


(9) 


^2 • \ Sin (u)t + .j) - -£^) (10) 

63 - E„ sin (a)t + ♦ + Jl-) ( 11 ) 

If equations ( 6 ) - (11) are substituted Into (5), the simplified 
result for Instantaneous developed motor torque Is found. 

Td - 7 EmI I I- cos (H.) {cosC(n-l)J mt - (^3 - co8C(n+l)i a.t + ♦3X12) 
n “ o o 6 

n - 1,5,7,11,13,*** 


From (12), It Is seen that the instantaneous developed motor torque 
Is made up of a constant term plus all multiples of the sixth harmonic 
of motor phase current frequency. Further, the multiples of sixth 
harmonic components of torque are a direct result of the current harmonics 
that lie Immediately on either side of multiples of six. Specifically, 
the sixth harmonic of torque results from the fifth and seventh harmonic 
of current; and, the twelfth harmonic of torque results from the eleventh 
and thirteenth harmonic of current; etc. Clearly, the motor developed 
torque will be constant If the harmonics of motor phase current described 
by 6 n ± 1 , n " 1,2,3,••• are eliminated. 

The above dlsctissed harmonics only need to be eliminated If their 
existence results In torque harmonics in the range of response of a 
coupled mechanical load. Since typical mechanical loads esdiiblt negligible 
response to frequencies above 20 hertz. It should only be necessary to 
eliminate those nontrlplen, odd current harmonics beyond the first for 
which 


M < 


80 TT 

6 P“ni 


(13) 


where p number of poles and Is the motor speed In rad/s. 

5.0 MODULATION OF PHASE VOLTAGE 

The cycloconverter system Is voltage excited; thus, if the (n ± 1) 
harmonics are to be eliminated from motor phase currents, they must be 
eliminated from the phase voltage set. 

5.1 Development of Modulation Function 

Modulation function h(ut) that contains no selected harmonics up 
through M Is derived In Appendix B. The balanced three-phase voltage set 
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to be applied for equations (1) - (3) Is given by 


^an 

- v^ h((»)t) 

(14) 

^bn 

2 it 

- Vj h(u»t - ) 

(15) 

^cn 

- v^ h(a)t + ^ ) 

(16) 


If Is the radian frequency of the three-phase AC source, then Is a 
periodic wave form described by 


Vd(t) - 


and, Vd(t + T) ■ 


where T ■ 


It Is shown In Appendix A that If the modulated phase voltages 
('''an»^nt''’'cn) fo contain none of the selected harmonics through M, 

It Is necessary that the freqxiency of v^ be much greater than the 
frequency of h(t); however, this condition Is easily satisfied at near 
zero speed. 

Figure 4 llliistrates the unmodulated phase voltage. The modulated 
phase voltage Is shown by Figure 5. 

5.2 Minimization of Slack Variable 

Development of the modulation function for use with the synchronous 
envelope cycloconverter In Appendix A required an extension of the well- 
known pulse width modulation techniques for elimination of harmonics In 
Inverter drives Cl4] to Include a slack variable. The modulation function 
contained no selected harmonics through M as determined by (13). For 
the reported work, the slack variable Is always chosen so as to minimize 
the sum of the squares of the amplitude of the next higher ( 6 n ± 1 ) 
current harmonic pair beyond harmonics selectively eliminated. Listings 
of the FORTRAN programs used to compute the modulation function are 
presented In Appendix B. 


slnCw^t + g- + «s) » 0 < Wj^t < — (17) 

vj(t) ( 18 ) 

iL. 
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5.3 Cycloconverter Connection Method 


Since the modulation ftmctlon hCut) does not repeat Itself every ir/6 
radians on the Interval Cir/6, Sir/63, attention must be directed to the 
cycloconverter connection scheme utilized so that the eliminated harmonics 
are not Inadvertently reintroduced. Specifically, a conduction path must 
be provided so that all harmonics remaining In the phase voltage can flow 
If the phase current Is to be of the same harmonic content. Two connection 
realizations are possible wherein no restriction on phase current harmonics 
are Imposed by Klrchhoff's current law. The first of these acceptable 
connections Is that of the midpoint cycloconverter with the neutral lead 
In place. A second realization would be the full-bridge cycloconverter 
with total phase Isolation. Since this latter connection requires 36 SCRs, 
It was not used In the reported study due to a factor of three on component 
requirement over the midpoint cycloconverter case. However, It should be 
pointed out that the full-bridge cycloconverter with Its six-pulse output 
has a smaller ripple magnitude and may offer a slight efficiency advantage. 

5.4 Fourier Spectmm Analysis 

Nature of the waveforms for current and Instantaneotis developed 
torqiie are such that visual Inspection to determine harmonic content Is 
not practical or accurate. The listing of a FORTRAN program to find the 
Fourier coefficients Is given In Appendix D. 

This program was used to generate the data for the normalized Fourier 
spectrum plots of Section III for harmonics up through forty-eight In order 
to verify thap attempt to selectively eliminate harmonics was successful. 




















SECTION III 
RESULTS 


1.0 SUMMARY 

Three speed condltloas were selected for nunerlcal study to evaluate 
the effectiveness of the harmonic elimination technique In removing 
undeslred torque pulsations from the midpoint cycloconverter driven 
brushless DC motor system when using synchronous envelope control: 

(a) . 4500 rpm (lOZ speed case) 

(b) . 2250 rpm ( 5Z speed case) 

(c) . 450 rpm ( IZ speed case) 

Although these speeds are great enough that elimination of harmonics as 
determined by equation (13) would not be necessary, there are two reasons 
for their selection. The first Is conservation of computer time. Computa¬ 
tion of a set of data for the 450 rpm (IZ speed case) required approxi¬ 
mately 8 hours of CPU time on a Control Data PDP 11/44 computer. Since 
many trlal-and-error runs were necessary for each data set, runs at lower 
motor speed become prc^bltlve timewise. Secondly, If the harmonics can 
be successfully ellMnated at higher speeds, then there Is not difficulty 
to be encountered at reduced speed. Justification lies In the fact that 
the undeslred harmonics are totally eliminated from the phase voltage set. 
Any re-lntroductlon of the undeslred harmonics Into the phase current 
wave forms would be attributable to commitation overlap which is dependent 
on the La/R time constant that determines current decay at the end of 
each phase current conduction sequence. Since this time constant Is not 
speed dependent, the conmutatlon overlap interval Is a larger portion of 
the current wave form period at high speed than It Is at low speed. Thus, 
the departure of phase current from the form of the phase voltage Is 
greater at high speed than at low speed. It Is concluded that If 
undeslreable harmonics are not re-lntroduced at high speed, then success 
at low speed Is guaranteed. 

For each speed case, the nature of current and torque were examined 
for application of a balanced set of phase voltages of the following 
nature: 

(a) . No harmonic eliminated. 

(b) . Fifth and seventh harmonics eliminated. 

(c) . Fifth, seventh, eleventh, and thirteenth harmonics eliminated. 
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A brushless DC machine system has an extra Input or degree of freedom 
over either a conventional DC machine or a s 3 nachronous machine drive 
system In that the angular displacement between the stator and rotor 
magnetic fields can be specifically controlled. For this study, this 
extra degree of control freedom Is Implemented by specification of the 
angle (a) measured from the onset of positive phase current conduction 
to the positive going zero crossing of the associated motor phase counter 
emf. A convenient value (a - 30^) Is used for most of the study; however, 
since selection of a shifts the position of motor phase counter emf with 
respect to phase current, the nature of coimnutatlon overlap can be altered 
by the value of a. The effect of changing a over the range from 15° to 
40° Is studied for one speed case. 

2.0 PERFOBMANCE STUDIES 


Eqtiatlons (1) - (3) can be written In compact matrix form as 


^1 . ^CR3l.^^-^e 


( 20 ) 



"ir 

(21) ; e - ^ 
” P 

sln(u)t-oX 

1. - 

^2 

sln(ut-o- ) 

3 




_sln(a)t-a+ ^ ) 


( 22 ) 


CR3 - 

Ra 

(23) ; v^ - v^(t) 

h(u)t) 

h(u.t - ^ ) 


R 

L 3 J 


h(u)t + ) 


(24) 


Equation (20) Is a nonlinear set of differential equations due to the 
elements of CR] being functions of current. Any attempt at linearization 
for a wide range of ^ would unacceptably alter the physical nature of 
the problem by removing the rectifying characteristic of the SCRs. Conse¬ 
quently, a closed form solution of (20) Is not possible and a numerical 
Integration must be Implemented. A fixed Increment (1.0 x 10~^s) fourth- 
order Runge-Kutta numerical Integration method Is utilized for this work. 

After forming a modulation function through use of the programs of 
Appendix B, selecting a particular motor speed (Ug^ > 2u/p), and arbitrarily 
choosing a valtie of SCR firing angle Og(see equation (17)), numerical 
solution to (20) Is Implemented subject to Initial conditions ^(0) > 0^. 
Integration Is continued until steady-state Is reached. Instantaneotis 
values of developed torque are calculated according to (5) and numerically 









averaged over the last steady-state cycle of integration to give the 
average value of developed torque. This value of average torque is 
compared with a specified value for the operating point. If the average 
value calculated is different from the specified value, then Ug is 
appropriately adjusted and the numerical integration is repeated until the 
difference magnitude between calculated and specified torque are suffi¬ 
ciently close in an epsilon sense. For all cases studied, the specified 
torque was taken as 1.5 N*m. 

The listing of a FORTRAN program to compute the performance as 
described above is presented in Appendix D. In addition to the discussed 
computations, the program also calculates the RMS value of phase current 

over the last steady-state cycle as ( ^ ^ i^ At)^. Also, values of phase 

current and Instantaneous torque are stored over the last steady-state 
cycle for use by the Fourier Spectrum Program of Appendix C. 

Results of the various speed cases and harmonic elimination studies 
are summarized in the paragraphs that follow. 

2.1 4500 RFM (lOZ Speed Case) 

2.1.1 Unmodulated Phase Voltage . 

Input Data 

Motor speed, n^j « 4500 rpm 
Motor frequency, f ■ 150 Hz 
SCR firing angle, Og - 58.68° 

Specified average torque, Tjg^ - 1.5 N*m 


Output Data 

Average torque, Tjav “ 1.501 N*m 
RMS phase current, Ijms ■ 72.77 A 


A plot of phase current ij is given by Figure 6. The associated 
Fourier spectrum of 1^^ is shown by Figure 7. Predominantly trlplen 
harmonic current that flows through the neutral connection is depicted by 
Figure 8. 

A plot of the instantaneous developed torque where the sixth harmonic 
is obvious is given by Figure 9. A Fourier spectrum of the instantaneous 
torqtie is shown in Figure 10 where existence of all multiples of six as 
theoretically predicted are noted. 
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Figure 7. Fourier Spectrum of Phase Current - A500 RPM, No Modulation 
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Figure 10. Fourier Spectrum of Developed Torque - 4500 RPM, No Modulation 







2.1.2 Elimination of Sixth Harmonic of Torque . The sixth harmonic of 
torque Is eliminated by removing the fifth and seventh harmonics of phase 
current. 


Motor speed, n^ ■ 4500 rpm 

Motor freqtiency, f ■ 150 Hz 

SCR firing angle, - 46.45° 

Specified average torque, - 1.5 N*m 
Angles of modulation function, a. - 40.76° 
oi - 47.73° 
03 - 58.65° 


Average torque, T, « 1.5001 N*m 

RMS phase current, I^j^g - 74.24 A 


A plot of phase current Ij^ Is given by Figure 11. The associated 
Fourier spectrum of 1^ Is shown by Figure 12, where negligibly small fifth 
and seventh harmonics are noted. Current flow through the neutral line 
la depicted by Figure 13. 

A plot of the Instantaneous developed torque Is given by Figure 14. 
The Fourlet spectrum of the developed torque Is shown In Figure 15 where 
It Is seen that the sixth harmonic Is near zero. 

2.1.3 Elimination of Sixth and Twelfth Harmonics of Torque . Both the 
sixth and twelfth harmonics of torque are eliminated by removing the fifth, 
seventh, eleventh, and thirteenth harmonics of phase current. 


Motor speed, n^^ « 4500 rpm 

Motor frequency, f ■ 150 Hz 
SCR firing angle, a, - 42.63° 

Specified average torque, Tj^y • 1.5 N>m 
Angles of modulation function, ai ■ 38.73° 
02 - 42.13° 
03 - 57.25° 
04 - 61.93° 

05 - 66.86° 


Average torque, T^av " 1*505 N.m 

RMS phase current, Ig^S ~ 74.59 A 
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Figure 12. Fourier Spectrum of Phase Current - 4500 RPM, Fifth and 
Seventh Harmonics Eliminated 
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A500 RPM, Fifth and Seventh Harmonics Eliminated 
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A plot of phase current Is given by Figure 16. The associated 
Fourier spectrum of 1^ Is shown by Figure 17 where negligibly small fifth, 
seventh, eleventh, and thirteenth harmonics are noted. Current flow 
through the neutral line Is depicted by Figure 18. 

A plot of the Instantaneous developed torque Is given by Figure 19. 
The Fourier spectrum of the developed torque Is shown In Figure 20 where 
It Is seen that the sixth and twelfth harmonics are near zero. 

2.2 2250 RFM (5Z Speed Case) 

2.2.1 Unmodulated Phase Voltage . 

Input Data 

Motor speed, n^ * 2250 rpm 
Motor frequency, f - 75 Hz 
SCR firing angle, 3 - 61.93° 

Specified average torque, 1 ^ 3 ^ > 1.5 N*m 

Output Data 

Average torque, Tjav “ 1*^99 N*m 

RMS phase current, Ipjjs “ 69.53 A 


A plot of phase current l^ Is given by Figure 21. The associated 
Fourier spectrum of 1^ is shown by Figure 22. 

A plot of the Instantaneous developed torque Is given by Figure 23. 

A Fourier spectrum of the developed torque Is shown by Figure 2A. 

2.2.2 Elimination of Sixth Harmonic of Torque . The sixth harmonic of 
torque Is eliminated by removing the fifth and seventh harmonics of phase 
current. 


Input Data 

Motor speed, nm * 2250 rpm 

Motor frequency, f - 75 Hz 
SCR firing angle, Og - 50.21° 

Specified average torque, T^gv * 

Angles of modulation function, a. - 40.76*^ 
02 - 47.73° 
03 - 58.65° 


Output Data 

Average torqiie, ■ 1.5001 N-m 

RMS phase current, IpMS ~ 73.96 A 
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Figure 16. Phase Current - ASOO RPM, Fifth, Seventh, Eleventh 
and Thirteenth Harmonics Eliminated 













Figure 17. Fourier Spectrum of Phase Current - 4500 RPM, Fifth, 

Seventh, Eleventh, and Thirteenth Harmonics Eliminated 
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A plot of phase current Ij^ is given by Figure 25. The associated 
Fourier spectrum Is shown by Figure 26. 

A plot of the Instantaneous developed torque Is given by Figure 27. 

The Fourier spectrum of the developed torque is shown by Figure 28. 

2.2.3 Elimination of Sixth and Twelfth Harmonics Of Torque . Both the 
sixth and twelfth harmonics of torque are eliminated by removing the fifth, 
seventh, eleventh, and thirteenth harmonics of phase current. 

Input Data 

Motor speed, Ujjj - 2250 rpm 

Motor frequency, f ■ 75 Hz 
SCR firing angle, Og * 47.31° 

Specified average torque, T^j^^ » 1.5 N*m 
Angles of modulation function, ai ■ 38.73° 

02 - 42.13° 

03 - 52.25° 

04 - 61.93° 

05 - 66.86° 

Output Data 

Average torque, T^j^^ » 1.5003 N*m 
RMS phase current, l{0{s * 74.64 A 

A plot of phase current l^ Is given by Figure 29. The associated 
Fourier spectrum of Ij^ Is shown by Figure 30. 

A plot of the Instantaneous developed torque Is given by Figure 31. 

The Fourier spectrum of the developed torque is shown by Figure 32. 

2.2.4 Elimination of Sixth and Twelfth Harmonics of Torque (a " 15°) . The 
angle from the onset of positive phase current to the positive going 
crossing of the associated motor phase coxmter emf Is changed from 30° to 
15° to examine effect on current harmonics. 


Input Data 


Motor speed, ng, * 2250 rpm 
Motor frequency, f • 75 Hz 
SCR firing angle, a, > 59.84° 
Specified average torque, T^gy • 
Angles of modulation function, ai 
02 

03 

“4 

“5 


1.5 N*m 

- 38.73° 

- 42.13° 

- 52.25° 

- 61.93° 

- 38.73° 
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Figure 26. Fourier Spectrum of Phase Current - 2250 RPM, Fifth 
and Seventh Harmonics Eliminated 
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Fourier Spectrum of Developed Torque - 2250 RPM, 
Sixth Harmonic Eliminated 
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Fourier Spectrum of Developed Torque - 2250 RPM, 
Sixth and Twelfth Harmonics Eliminated 





Average torque, ■ 1.502 N*m 

RMS phase current, - 53.33 A 


A plot of phase current 1. Is given by Figure 33. The associated 
Fourier spectrum is shown by Figure 34. 

A plot of Instantaneous developed torque Is given by Figure 35. 

The Fourier spectrum of the developed torqvie Is shown by Figure 36. 

The re-introduced harmonic magnitudes are reduced to a lower level 
than for a - 30°. 

2.2.5 Elimination of Sixth and Twelfth Harmonics of Torque (a ■ 40°) . 
The angle from onset of positive phase current to the positive going 
crossing of the associated motor phase counter emf Is set at 40° to 
examine the effect on current harmonics. 

Input Data 

Motor speed, Oq « 2250 rpm 
Motor frequency, f « 75 Hz 
SCR firing angle, Ug - 16.77° 

Specified average torque, Tjgy - 1.5 N-m 
Angles of modulation function, ■ 38.73° 

02 - 42.13° 

03 - 52.25° 
o^ - 61.93° 

05 - 66.86° 


Output Data 


Average torque, T^g^ “ 1*5003 N*m 
RMS phase current, IpMS ” 107.52 


A plot of phase current 1^^ is given- by Figure 37. The associated 
Fourier spectrum of phase current Is shown by Figure 38. 

A plot of Instantaneous developed torque Is given by Figure 39. 

The Fourier spectrum of developed torque is shown by Figure 40. 

Although the angle power factor angle (o i- 30°) has been extended to 
75° resulting In a large current to maintain the same power flow, no 
significant re-lntroductlon of undeslred harmonics Is noted. 
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Phase Current - 2250 RPM, a - IS®, Fifth. Seventh, 
Eleventh, and Thirteenth Harmonica Eliminated 
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Figure 37. Phase Current - 2250 RPM, a - AO®, Fifth, Seventh, 
Eleventh, and Thirteenth Harmonics Eliminated 
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2.3 450 RFM (IZ Speed Case) 

2.3.1 Unmodulated Phase Voltage . 

Input Data 

Motor speed, nm > 450 rpm 

Motor freqtiency, f ■ 15 Hz 
SCR firing angle, Og - 64.52° 

Specified average torque, N'® 

Output Data 

Average torqvie, T^j * 1.4996 N*m 
RMS phase current, Ipjjs “ ^ 


A plot of phase current li Is given by Figure 41. The associated 
Fourier spectrum is shown by Figure 42. 

A plot of Instantaneous developed torque is given by Figure 43. A 
Fourier spectrum of developed torque is shown by Figure 44. 

The .instantaneous wave form plots for the low speed case do not dis¬ 
play all of the rippled due to the lack of capacity of the plot routine 
to store sufficient points. 

2.3.2 Elimination of Sixth Harmonic of Torque . The sixth harmonic of 
torque is eliminated by removing the fifth and seventh harmonics of phase 
current. 


Input Data 

Motor speed, Oq ~ 450 rpm 

Motor frequency, f * 15 Hz 

SCR firing angle, * 52.53° 

Specified average torque, - 1.5 N’m 

Angles of modulation function, a. - 40.76° 
02 - 47.73° 
03 - 58.65° 

Output Data 

Average torque, T^j^y * 1.4996 N*m 
RMS phase current, Irms “ 75.2 A 


A plot of phase current 1^ is given by Figure 45. The associated 
Fourier spectmm is shown by Figure 46. 
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A plot of Instantaneous developed torqxie Is given by Figure 47. 
The Fourier spectrum of developed torque Is shown by Figure 48. 

2.3.3 Elimination of Sixth and Twelfth Hantonlcs of Torque . Both the 
sixth and twelfth harmonics of torque are eliminated by removing the 
fifth, seventh, eleventh, and thirteenth harmonics of phase current. 

input Data 


Motor speed, n^ « 450 rpm 
Motor frequency, f » 15 Hz 
SCR firing angle, Og - 49.78°^ 
Specified average torque, T^^^ 
Angles of modulation function. 


- 1.5 N*m 


“1 

02 

“3 

“4 

“5 


38.730 

42.13° 

52.25° 

61.93° 

66 . 86 ° 


Output Data 

Average torque, T^^ ■ 1.507 N»m 
RMS phase current ipjjs ■ 77.38 A 


A plot of phase current Is given by Figure 49. The associated 
Fourier spectrum of phase current Is shown by Figure 50. 




A plot of Instantaneous developed torque Is given by Figure 51. The 
Fourier spectrum of developed torque Is shown by Figure 52. 
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Figure 46, 

















































100 



Harnonic Nunber 

Figure 50. Fourier Spectrum of Phase Current - 450 RPM, Fifth, 

Seventh, Eleventh, and Thirteenth Harmonics Eliminated 
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Developed Torque - 450 RPM, Sixth 
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lectrum of Developed Torque - 450 RPM, 
Twelfth Harmonics Eliminated 










SECTION IV 


DISCUSSION 


1.0 INCREASE IN OHMIC LOSSES 

Since modulation of the wave form definitely redistributes the 
harmonic content of the phase current wave form, an examination of the 
change In winding ohmic losses Is In order. For the 4500 rpm case It 
Is noted that the RMS phase current Is found to be as follows: 

• ^BMS * 72.77 A for no modulation 

^RMS " 74.24 A for fifth and seventh harmonic eliminated 

(c). Ip^ o 74.59 A for fifth, seventh, eleventh, and thirteenth 
harmonic eliminated 

The Increase In ohmic losses assuming a constant winding resistance Is 
given by the ratio of RMS current for a changed condition to the RMS 
value of current for the unmodulated wave form. It Is concluded that 
the ohmic losses are Increased by 2.1Z for the case of elimination of 
fifth and seventh harmonics. The ohmic losses are Increased by 2.5Z for 
the case of elimination of the fifth, seventh, eleventh, and thirteenth 
harmonics. 

Calculated data for the 2250 rpm and 450 rpm cases show Increases 
of 7 to 15Z; however, there Is no theoretical basis for greater Increase 
In these cases. The reason for this apparent larger Increase for the 
lower speed cases Is a numerical problem due to storing and calculating 
the RMS valtie based on a limited number of data points (2000) for trans¬ 
mittal to the plot routine. As a result, the high frequency ripple 
currents are sampled at Irregular Intervals giving rise to Inaccuracy In 
calculation of the RMS value. For the case of 4500 rpm, the data point 
limit allows storage of several points during each ripple cycle, thus, the 
resulting RMS values are considerably more accurate than for the other 
two speed cases. 

2.0 CYCLOCONVERTER MODE CHANGE 

As discussed and Illustrated by Figure 2, a neutral connection Is 
necessary to successfully eliminate selected phase current harmonics. 

Since harmonic elimination Is only necessary at speeds near zero, and 
since an ohmic loss penalty Is Incurred as a result. It might be 
deslreable to not operate with FWM control over the full range of speed. 
The Harmonic Eliminator Switch In the neutral line of Figure 2 could be 
opened for higher speeds and the cycloconverter mode of operation changed 
to phase-control from synchronous envelope. 
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3.0 mENSION TO HIGHER HABMONIC ELIMINATION 

For specific Illustration, this study has only considered elimination 
of torque harmonics through the twelfth. As motor speed approaches zero, 
it may be necessary to eliminate even higher torque harmonics. Theoreti¬ 
cally, there is no limitation to be encountered. 
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SECTION V 


CONCLUSIONS AND RECOMMENDATIONS 


The work performed by this study shows that systematic reduction of 
harmonic torque pulsations In brushless DC drives Is feasible by use of 
selective current harmonic elimination. Further, there appears to be 
only a small ohmic loss penalty. 

Further theoretical study should be made to properly assess quanti¬ 
tatively the effect of Increased motor leakage Inductance and of Interphase 
reactors on "reintroducing" the harmonics eliminated from the phase 
voltage due to extending the commutation overlap. 

Also, a hardware realization of a midpoint cycloconverter brushless 
DC drive using the control principles formulated In this study should be 
made to verify the results. 
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APPENDIX A 


HARMONICS OF PWM TIME FUNCTION 


A modulation function h(ut) was derived In Appendix B and used In 
Section II - 5.0 to form a modulated phase voltage. For example. 

Van " h<“t) vjCt). 

Harmonic content of a modulated voltage Is of concern. The study 
Is best understood by use of a typical example. Let h(ut) contain all 
odd harmonics except the fifth and seventh, then 


h(ut) « a^cos wt + a^cos 3<Dt + a^cos 9«t + a^jcos Hut + ••• 

+ a„cos nut + ••• (A.l) 


Assume the Vj(t) contains a DC term plus one high frequency harmonic, then 

Vj(t) ■ bg + bjjjCos mut (A.2) 

Forming v^ as the product of h(ut) and Vj(t) and slnqillfylng yields 

Vj(t) ■ c^cos ut + C3COS 3«t + C3C0S 9ut + c^^cos Hut + ••• 

+ dj^j^cos(m+l)ut + d^_j^cos(m-l)ut + d^2Cos(m+3)ut 
+ dj,_3COs(m-3)ut + ••• (A. 3) 


It Is seen from (A.3) that through the cross product terms, the PWM 
expression may contain frequencies that were selectively eliminated from 
h(ut). However, the coefficients (dg^i, dB_i, •••) of (A.3) are products 
of Sq and bn. Since Sq decreases as 1/n, If m Is large, then the 
coefficients of the low frequency components of h(ut) resulting cross 
product terms will be negligibly small. It Is concluded that the fre¬ 
quency of the AC source must be large compared with the motor frequency 
If the PWM function Is to not contain the harmonics selectively eliminated 
from the modulating function h(ut). But, at near zero speed, the source 
frequency Is large compared with the motor frequency. 









APPENDIX B 



r:l 


MODULATION FUNCTION PROGRAMS 


1.0 INITIAL SOLUTION 

As first procedure, a candidate or Initial modulation function, h(u>t), 
is found to serve as a basis for an optimization search to determine the 
final modulation function. 

1.1 Theory of Harmonic Elimination 

Extending Patel and Hoft's pulse width modulation technique Cl4D to 
Include a slack variable, it is possible to find a modulation function 
that has selected harmonics elimination while existing on the quarter wave 
interval from 30° to 90°. 

The Fourier series representation of the modulation function assuming 
odd quarter wave symmetry is 

f((Dt) * \ a^slnCnut) , for odd n 

where 

®M 

a„ ■ /qj sln(nut)dut + sln(nut)dut + ••• sln(nu>t)dut, (B.l) 


30° < Oj <02 < • • • < ojj < 90° 


This reduces to 


I 


^ (-l)^+^cos(na^) (B.2) 

for M both odd and even. The number of harmonics eliminated is in general 
the same as the number of switching transitions per quarter cycle. For M 
greater than 2 the equations cannot be solved in closed form, so numerical 
methods are needed to obtain solutions. The following normalized system 
of equations Is to be used to obtain solutions: 


^nl " I (-•l)*''’^cos(nia|j) - 0 

^n2 " I (-l)‘"‘^^cos(n2a^) - 0 


^nM “ I (-l)*^'‘'^cos(njja|j) » 0 










where ni,n2>***i^ are the harmonics to be eliminated. 


Many of the solutions of (B.3) are meaningless or Impractical, so it 
is necessary to search by brute force using a loose convergence criterion 
for candidate solutions, and then to use a numerical method such as Newton 
Raphson to obtain a refined solution. 

The basic method of selecting candidate modulation function for 
further optimization Is Illustrated by the flow chart of Figure B.l. A 
candidate set Is one that meets the following criteria: 


^squar ^ ^min 

^squar " ^nl ^n2 *“ ^nM 

If for the candidate modulation function, 30® < a^ < a2 < * * * < aj^ < 90°, 
then the solution Is considered for optimization. Otherwise, the solution 
Is discarded. 








SELECT M HARMONICS 
AND SLACK HARMONIC 
TO BE ELIMINATED 



Figure B.l. Flow Chart for Selection of Candidate 
Modulation Functions 












1.2 Initial Solution Program 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c c 
c main program c 


ccecceeccccccccecccccccccccccccccccccceccccccccccccccccccccccccccccccccc 
ccccccccccccccccccoccccccccccccccccccccccccccccccccccccccccccccccccccccc 
on- number of harmonics to be eliminated 
c this includes the slack variable 

- this is the amount that each alpha 
will be incremented by during the 
infinite search 

- newton rapson method is used to 
search for a solution if the sum 
of squares is less than fmin 

> these are the harmonics to be reduced 

- this is the minimum value of alpha 
for the infinite search 

- not used 

- not used 

dimension a(10},r(10),f(10) 
coonon pl,conv,lwr,del 
real inc 
call sopen(O.O) 
pla3.141592653 


inc 


fmin 


r(l) 

amin 


iwr 

del 


kIminaO 

k2min30 

k3mins0 

kUmlnsO 

k5min30 

kbminsO 

kl3i 


k23l 

k3s1 

k4sl 

k53l 

k63l 

conv3 180. 0/pi 
read(9»*) n 
write(6,») ' ns»,n 
read(9»*) inc 
write(6,») ’ ine3',inc 

read(9t*) flain 

wrlte(6,») ’ fmln3*,fmln 

read(9.») (r(l),l3i,7) 

write(6,*) ' harmonics3',(r(i),i3l,n) 

read(9t*) amin 

write(6,*) ' amins'^amln 

read(9»*) Iwr 

wrlte(6,*) ' iwr3»,iwr 











maxslnt( (90.0<-afflln)/inc ) 

k1maxsint( (90.0-amln)/lnc ) 

k2maxsmax 

kSmaxsmax 

kUmaxsmax 

kSmaxsmax 

k6maxsmax 

IsepsintC 1.0/inc ) 

cccccccccccccccccccccccccccccccccccceccccccccccccccccccccccccccccccccccc 

e 

c tha following are essentially do loops but were 

c necessary to get past a problem with the compiler 

c 

kimlnsi 

100 k1sk1mln-1 

105 klakl+1 

lf( n .eq. 1 ) go to 1000 

k2mlnsk1-»>lsep 
200 k2ak2mln-1 

205 k2sk2-i'1 

lf( n .eq. 2 ) go to 1000 

k3mlnsk24-lsep 
300 k3«k3nln-1 

305 k3«lc3>1 

lf( n .eq. 3 ) go to 1000 

k4mlnsk3’*’lsep 
400 k4sk4mln-1 

405 k4sk4<i>1 

lf( n .eq. 4 ) go to 1000 

k5mlnxk4-i>lsep 
500 k58k5mln-1 

505 k5sk5+1 

lf( n .eq. 5 ) go to 1000 

k6mlnsk5-*>l8ep 
600 k6sk6mln-1 

606 k6sk6-f1 

lf( n .eq. 6 ) go to 1000 

1000 a(1)*( float(k1}*lnc4-amln )/conv 

a(2)«( float(k2)*lnc-^amln )/oonv 
a(3)»( float(k3)*lno4>amln )/oonv 
a(4)s( float(k4)*lnc-«-amln )/conv 
a(5)s( float(k5)*lno-»>amln )/oonv 
a(6)s( float(k6)*lnc4>amln )/eonv 
call ff( a,r,f,n ) 
call square( ftfsquar.n ) 

lf( faquar .It. fmln } call look( a,r,f,fsquar,n ) 
lf( k6 .It. k6max .and. n .ge. 6 ) go to 606 

lf( k5 .It. k5max .and. n .ge. 5 ) go to 505 

lf( k4 .It. k4max .and. n .ge. 4 ) go to 405 

lf( k3 .It. k3niax .and. n .ge. 3 ) go to 305 
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3000 

3100 

3200 


if( k2 .It. k2max .and. n .ge. 2 ) go to 205 
if( k1 .It. kimax .and. n .ge. 1 ) go to 105 
wrlte(6,3100) 

wrlte(6,3000) ( a(i)*conv,is1,n ) 

call sclose(O.O) 

formate 10(e15.7.5x) ) 

formate* as*) 

formate* harmonicas') 

stop 

end 
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ccoccccccoccccccccccccccccccceceeccccccccccccccccccccccccccccccccccccccc 
coccccccccccccocccccccccccccceccececcccccceccccccccccccccccccccccccccccc 
c c 

c subroutine a c 

c c 

cccccccceeccocceccccoccoccccccccccceeccccccccceccecccccccccccccccecceecc 
cccccccccooeococcccececcccoccccccccccccccccccccccccccccccccccccccccccccc 
subroutine fa( a,da,n ) 
difflension a(10),da(10) 
common pltConv,lwr,del 

do 200 Isl.n 
a(i)sa(i)-da(l) 

200 continue 
return 
end 








cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c c 
c subroutine close c 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
subroutine sclose(duimny) 
closed) 
close(2) 

closeO) 

close(4) 
close(7) 
close(8) 
close(9} 
return 













cccccccccccccccccccccccccccccccccccccccccccccocccccccccccccccccccccccccc 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c c 

c subroutine f c 

c c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
subroutine ff( a,r,f,n ) 
dimension a(10),r(10),f(10) 
common pl,conv,lwr,del 
do 200 ls1,n 
signsl.O 
f(i)=0.0 

do 100 Jsl.n 

f(i)sf(i)+algn«coa( r(i)»a(j) ) 
signs-sign 

100 continue 

200 continue 
return 
end 
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cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


c c 

c subroutine look c 

c c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
ccccecccccccccoccccceccccccccccccccccccccccceocccccccccccccccccccccccccc 
subroutine look( a,r,f,fsquartn ) 

dimension da(10),r(10),a(10),f(10),df(10,10),dfi(10,10) 

common pi,conv,iwr,del 

write(6,3050) 

write(6,3050) 

write(6,3100) 

write(6,3000) ( a(i)*conv,i=1,n ) 
write(6,3150) 

write(6,3000) ( f(i),i=1,n ) 

write(6,3200) 

write(6,3000) fsquar 

kounts30 

fminsO.01 

delta=0.001 


do 2500 loopsl,kount 
call fdf( a,r,f,df,n ) 
call matv( df,n,det,dfi ) 
call multK f,dfi,da,n,n ) 
call fa( a,da,n ) 
call ff( a,r,f,n ) 
call square( f,fsquar,n ) 

lf( loop .gt. 5 .and. fsquar .It. fmin ) go to 2900 
2500 continue 

2900 write(6,3050) 

write(6,3100) 

write(6,3000) ( a(i)*conv,is1,n ) 
write(6,3150) 

write(6,3000) ( f(i),isl,n ) 
write(6,3200) 
write(6y300&) fsquar 
3000 formate 10(e10.4,5x) ) 

3050 formate » • ) 

3100 formate ' as* ) 

3150 formate ’ f=' ) 

3200 formate ' fiquars’ ) 

4000 formate'test',i3) 

5000 return 

end 
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ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


c c 

C SUBROUTINE MATV C 

c c 

c c 


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE MATV(A,N,DET,B) 

DIMENSION INDEX(20,2),IPVOT<20),A(10,10),B(10,10),PIVOT(20) 

ONEsl.O 

ZEROsO.O 

c 

DO 10 Isl.N 

DO 5 Jsl.N 
B(I,J)3A(I,J) 

5 CONTINUE 

10 CONTINUE 
DETsONE 
DO 15 Jsl.N 
15 IPVOT(J)sO 

DO 75 Is1,N 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


c c 

C FOLLOWING 15 STATEMENTS C 
C FOR SEARCH FOR PIVOT ELEMENTS C 
C C 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
TsZERO 
DO 40 Js1,N 

IF(IPV0T(J)-1) 20,40,20 
20 DO 35 Ks1,N 

IF(IPV0T(K)-1) 25,35,100 
25 CONTINUE 

DlsABS(T) 

D2sABS( B(J,K) ) 

IF(D1-D2)30,35,35 
30 IROWaJ 

ICOLsK 
TsB(J,K) 

35 CONTINUE 

40 CONTINUE 
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IPV0T(IC0L)=IPV0T(IC0L)+1 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


c c 

C FOLLOWING 10 STATEMENTS PUT C 
C PIVOT ELEMENT ON DIAGONAL C 
C C 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

IF(IROW-ICOL) 45,55,45 
45 DETs-DET 

DO 50 La1,N 
TsB(IROW,L) 

B(IR0W,L)sB(IC0L,L) 

50 B(IC0L,L)sT 

55 INDEX(I,1)sIR0W 

INDEX(I,2)3lC0L 
PIVOT(I)=B(ICOL,ICOL) 

DETsDET«PIVOT(I) 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


c c 

C FOLLOWING 3 STATEMENTS TO DIVIDE PIVOT C 

C ROW BY PIVOT ELEMENT C 

C C 


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

B(ICOL,ICOL)sONE 
DO 60 L3l,N 

60 B(ICOL,L)aB(ICOL,L)/PIVOT(I) 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC''JCCCCCCCCCCCCCCCCCCCCCCCC 


c c 

C FOLLOWING 7 STATEMENTS TO REDUCE C 
C NON-PIVOT ROWS C 
C C 


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

DO 75 LIs1,N 
IF(LI-ICOL) 65,75,65 
65 T*B(LI,ICOL) 

B(LI,ICOL)sZERO 
DO 70 Ls1,N 

70 B(LI,L}sB(LI,L)-B(ICOL,L)n 

75 CONTINUE 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


c c 

C FOLLOWING 11 STATEMENTS TO C 
C INTERCHANGE COLUMNS C 
C C 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
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DO 95 I»1,N 
LsN>I-i>1 

IF(IMDEX(L,1)-IMDEX(L.2) ) 85,95.85 
JROWsINDEX(L,1) 

JC0LslNDEX(L,2) 

DO 90 Xal.N 

TsB(K,JR0W) 

B(K,JR0W)sB(K,JC0L) 

B(K,JC0L)sT 

CONTINUE 

CONTINUE 

RETURN 

END 









cccccccceccccccccccccccececccccccecccccccccccccccccccccccccccccccccccccc 

cocccccccccecccccccccccccccccccccccccccccccecccccceccecccccccceccccccccc 


c subroutine multiply c 

c c 

ccccccccceeecccccccccecccoceccececcccccceccecccccccccccceccccecccccccccc 
cceccccccccccccccecccecccocccccccccccccccccccccccccccccccccccccccccccccc 
subroutine multi( a,b«c,10,lc0 ) 
dimension a(10),b(10,10),c(10) 
do 200 Is1.10 
c(l)s0.0 
do 100 ksl.kO 

c(l)sc(l)-t>b(l.k)«a(k) 

100 continue 

200 continue 
return 
end 
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ccceccceeccccecccccccccccccceeeccecccccccccccccccccccccccccccccccccccccc 
coecccccccccccccccccecccceeccccccccccccceecccccccccccccccccccccccceccccc 
c c 

c subroutine open c 

e c 

cceccccccccccccccceccccccccccccccccccccecccccccccccccccccccccccecccceccc 
ccccccccccceeecccccccooccccccecccccccccccccccccccccceccccccccccccccccccc 
subroutine sopen(dunny) 
open (unit s9. files •* input") 
rewind(9) 
return 
end 














ceecccccccccccceeeceeecccecccccccccccccccccccccccccccccccccccccccccccccc 
eeceeoeceeeeeococcccoecccccccccccccccccccccccccccccccccccccccccccccccccc 
c c 

c subroutine square c 

c c 

cccoccececccccoeeeeeeeceecccceeeeecccccccceeceecceecccccceeececccccccccc 
cccececeeeeeecoeoooooceecoeccoceeccccccceccecccccceccccccccccecccccecccc 
subroutine square( f,fsquar,n ) 
dimension f(10) 
fsquarsO.O 
do 200 ls1,n 

fsquarsfsquar4-f(l)*f(l) 

200 continue 
return 
end 
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2.0 OPTIMIZATION OF MODULATION FUNCTION 
2.1 Discussion of Procedure 

A candidate solution for the modulation function Is selected from 
the output of the Initial Solution Program for optimization. The slack 
harmonic of the Initial solution Is changed In small Increments and fed 
to a Newton-Raphson routine to find the new set of solutions. A con- 
tlnuoiis set of o's results from this procedure having 'Vilce*' properties. 
First, the new set of solutions is well behaved In that meaningless sets 
of solutions are not encountered and convergence Is quick. Second, the 
procedure gives a "visual" picture of the behavior of the solutions In 
relation to the slack variable. 

Since Infinite number of modulation functions exist, some method Is 
needed to select the "best" solution. The criteria used in this study 
was minimizing the sum of the squares of the magnitudes of the next set 
of current harmonics to be eliminated beyond those of present concern. 
For example. If the fifth and seventh harmonics were being eliminated 
from h( t), then the procedure Is to minimize the sum of the square of 
the magnitudes of the eleventh and thirteenth harmonics. 

A flow chart Illustrating optimization of the modulation function 
Is given by Figure B.2. 

Figure B.3 Illustrates the form of the modulation function for the 
case of elimination of fifth and seventh harmonics and for the case of 
elimination of fifth, seventh, eleventh, and thirteenth harmonics. 








L Figure B.2. Flow Chart for Optimizing 

Modulation Function 

















2.2 Modulation Function Optloilzatlon Program 


cccocccccccocccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

eccceceecceccececccccccccccccccceccccccccccccccccccccccceecccccceccccccc 


cccccecccccccccccccceoccccececececcccccccccccecceccccccccccccccccccccccc 
ccocecccccccccccccccccccccececcccccceccceccceccccccccccccccccccccccccccc 
c n - number of harmonics to be eliminated 

c n<t-1 harmonics are eliminated 

c r(i) - harmonics to be eliminated 

c a(i) - these are the angles from the infinite 

c search .l.e. this is the seed 

c rmax - upper limit of slack variable 

c rmin - lower limit of slack variable 

c points - number of solutions or points to be saved 

c msq - number of harmonics to be squared and sumed 

c tr(i) - harmonics to be summed and squared excluding tr(1) 

c tr(1) - any misc. harmonic that the user wishes to view 

c incre - used for increasing or decreasing the slack variable 

c incresl increment the slack variable from min to max 

c IncreaO decrement the slhck variable from max to min 

c hmag - this is the sum of sqaures 

dimension a(10),r(10),f(10),h(10),tr(10) 

common pi,conv 

call sopen(0.0} 

pls3.141592653 

convx180.0/pi 

read(9,*) n 

write(6,») ' na’.n 

read(9,») (r(i),ia1,5) 

wrlte(6,*) ' harmonicas*,(r(i)»is1»n) 

read(9,*) (a(i),isl,5) 

wrlte(6,*) ' alphass*,(a(i),is1,n) 

read(9t*) rmax 

wrlte(6,*) ' rmaxs',rmax 

read(9»*) rmin 

write(6,*) ’ rmins',rmin 

read(9»*) points 

wrlte(6,*) ' pointss*^points 

read(9»*) msq 

write(6,*) ' msqs',msq 

read(9,*) (tr(i),isi,5) 

wrlte(6,*) ' hsquar harmonicas*,(tr(i),is2,nisq-«>1) 

read(9t*) incre 
write(6,*) ' incres',incre 

rinos( rmax-rmin )/points 
maxrsintC rmax/rinc ) 
minrsintC rmin/rino ) 
irmaxsmaxrHiinr 








do 50 is1,5 
a(i)sa(i)/conv 
continue 

do 100 irsininr,maxr 

lf( incre .eq. 1 ) r(n+l)sfloat(ir)»rinc 
if( Incre .eq. 0 ) r(n-f1)sfloat( maxr+minr-ir )»rinc 
call ff( a,r,f,n>1 ) 
call square( ftfsquar,n-»-1 ) 
call lookC a,r,f,fsquar,n-f1 ) 
call share a,tr,h,n+1,maq+1 ) 
call shsqC h,hsquar,msq ) 
hmagssqrtC hsquar ) 
writed,*) ( a(i)*conv,is1,n<i>1) 

Hrite(2,1000) h(1) 
wrlte(3,1000) h(2) 
write(4,1000) h(3) 
wrlte(7,1000) r(n+1) 
write(8.1000) hoag 
continue 
call sclose(O.O) 
formate el0.4 ) 
formate 10ee15.7,5x) ) 
formate* a*’) 
formate* harmonics^') 
stop 
end 









cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

cccccccceccccccccccccccccccccccccccccecccccccccccccccccccccccccccccccccc 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
cccccccccccecccccceccocccccceccccecceccccccccccccccccccccccccccccccccccc 
subroutine fa( a,da,n ) 
dimension a(10),da(10) 
common pl.conv 

do 200 i=1,n 
a(l)=a(i)-da(i) 

200 continue 
return 
end 









cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
cecccccccccccccccccccccccccccccccccceccccccccccccccccccccccccccccccccccc 
c c 

c subroutine close c 

c c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
ccccccccccccccccccccccccccccecccccccccccccccccccccccccccecccccccccceceec 
subroutine sclose(duany) 
closed) 
close(2) 
close(3) 
close(4) 
close(7) 
close(8) 
close(9) 
return 
end 










cccccccccccccccccccccccccccccccccccccccccccccceccccccccccccccccccccccccc 
cccccccccccccceccccccccccccecccccccccccccccccccccccceccccccccccccccccccc 
c c 

c subroutine df c 

c c 

cccceceecccocccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
ccccccececcccccccccccccccccccccccccccccccccccccccccccccccccccccccccccecc 
subroutine fdf( a,r,f,df,n ) 
dimension a(10).r(10},f(10),df(10,10) 
conmon pi,conv 


do 200 is1,n 
Signs 1.0 

do 100 Jsi.n 

df(i,j)=-sign«r(i)»sin( r(i)*a(J) ) 
signs-sign 

100 continue 

200 continue 
return 
end 
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ccccccccscccccccccccccccccccececcccccccccccccecccccccccccccccccccccccccc 
cccccecccccccccccceccccccccccccccceceecccecccecceecccccccecccccccccccccc 
c c 

c subroutine harmonics c 

c c 

cccccccccccccccccccccccccccccceecccccccccccceccccccccccccccccccccccccccc 
cccccccccccccccccccccecccccccccccceccccccccccceccccccccccccccccccccccccc 
subroutine shar( a,tr,h,n,m ) 
dimension a(10),h(10),tr(10) 
common pl,conv 

do 200 ls1,m 
Signs 1.0 
h(l)s0.0 

do 100 jsl.n 

h(l)sh(l)>tJ».0/( tr(i)*pl )«slgn»cos( tr(l)»a(j) ) 
slgns-slgn 

100 continue 

200 continue 
return 
end 









ccccccccccccccccccccccecececcccceceeeeeececceeccccccccccccccccccccoccccc 

ceccccccccccccccccccccceccocececcceceeoceeeccccccccccccccccccccccccccccc 







ccececccccccceeccecccccceccccccceccccccccccccccccccccccccccccccccccccccc 

ceococccccecceeeccccecceeccccccceccccccccccccccccccceccccccccccccccccccc 


c c 

c subroutine look c 

0 c 

ecccoeccececccccceecccecccccceecceececcceccccccccccccccceccccccceccccccc 
cceeoecccccccccecccccccccccccceccecceccccccccecccccccccccccccccccccccccc 
subroutine look( a,r,f,fsquar,n ) 
dimension a(10),r(10),f(10),df(10,10),dfi(10.10) 
eoamK>n pi.conv 
kounts30 
fminsO.01 

do 2500 loopsl.kount 
call fdf( a,r,f,df.n ) 
call matv( df,n,det,dfi ) 
oall multK f,dfi,da,n,n ) 
call fa( a,da,n ) 
call ff( a,r,f,n ) 
call square( f,fsquar,n ) 

if( loop .gt. 5 .and. fsquar .It. fmin ) go to 2900 
2500 continue 

2900 dummysO.O 

3000 formate 10(e15.7,5x) ) 

3050 formate *0» ) 

3100 formate ' as* ) 

3150 formate ' fs' ) 

3200 formate ' fsquars' ) 

4000 formate'test',13) 

return 
end 














ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c c 


SUBROUTINE MATV 


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE HATV(A,N,DET,B) 

DIMENSION INDEX ( 20 , 2 ), IPVOT(20),A(10,10},B(10, 1 0),PIVOT(20) 

ONES 1.0 
ZEROsO.O 


c 

DO 10 1=1,N 

DO 5 Js1,N 
B(I,J)sA(I,J) 

5 CONTINUE 

10 CONTINUE 
DETsONE 
DO 15 Js1,N 
15 IPVOT(J)sO 

DO 75 Is1,N 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


c c 

C FOLLOWING 15 STATEMENTS C 
C FOR SEARCH FOR PIVOT ELEMENTS C 
C C 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
TsZERO 
DO 40 Js1,N 

IF(IPV0T(J)-1) 20,40,20 
20 DO 35 Ks1,N 

IF(IPV0T(K)-1) 25,35,100 
25 CONTINUE 

DlsABS(T) 

D23ABS( B(J,K) ) 

IF(D1-D2)30,35,35 
30 IROWsJ 

ICOLsK 
TsB(J,K) 

35 CONTINUE 

40 CONTINUE 


110 









IPVOT(ICOL)sIPVOT(ICOL) 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


C FOLLOWING 10 STATEMENTS PUT C 
C PIVOT ELEMENT ON DIAGONAL C 
C C 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
IF(IROW-ICOL) 45,55,45 
45 OETs-OET 

DO 50 La1,N 
T>B(IROW,L) 

B(IR0W,L)3B(IC0L,L) 

50 B(ICOL,L)sT 

55 INDEX(I,1)sIR0W 

INDEX(I,2)sIC0L 
PIVOT<I)*B(ICOL,ICOL) 

DETaDET*PIVOT(I) 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


c c 

C FOLLOWING 3 STATEMENTS TO DIVIDE PIVOT C 

C ROW BY PIVOT ELEMENT C 

C C 


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

B(ICOL,ICOL)sONE 
DO 60 La1,N 

60 B(ICOL,L)aB(ICOL,L)/PIVOT(I) 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


c c 

C FOLLOWING 7 STATEMENTS TO REDUCE C 
C NON-PIVOT ROWS C 
C C 


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

DO 75 LIa1,N 
IF(LI-ICOL) 65,75,65 
65 TsB(LI,ICOL} 

B(LI,ICOL)sZERO 
DO 70 Lsl.N 

70 B(LI,L)*B(LI,L)-B(ICOL,L)«T 

75 CONTINUE 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


c c 

C FOLLOWING 11 STATEMENTS TO C 
C INTERCHANGE COLUMNS C 
C C 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 







M 


DO 95 Isl.N 

LsN-l4-1 

IF(INDEX(L.1)-INDEXa.2) ) 85,95,85 
JR0WsIN0EX(L,1) 

JC0LsIN0EX(L,2) 

DO 90 Ks1,N 

TsB(K,JROW) 

B(K,JROW)sB(K,JCOL) 

B(K,JCOL)sT 

CONTINUE 

CONTINUE 

RETURN 

END 









ccccccccccccccccccccccccccccccccccccccccccccccceccceececcceccccccccccccc 

cccccccccccccccccccccecccccccceoccccccccceceeecccecceccccccccecccccccccc 


subroutine multiply 


cccccccccccccecccccccccccccccccccccccccccccccccccceecccccccceeeecccccccc 
cccccccccccccccccccccccccccccccccccccccccccccccccccecceecceccccccccccccc 
subroutine multK a,b,c,10,k0 ) 
dimension a(10),b(10,10),c(10) 
do 200 Is1,10 
c(l)s0.0 
do 100 ksl.kO 

c(l)sc(l)+b(l,k)»a(k) 

100 continue 

200 continue 
return 
end 
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eceecoccccoccccccccccccccccccccccccccccccccccccccccccccccccccccccecccccc 

cccececccccoccccccccccccccccccccccccoccccccccccccccccccccccccccccccccccc 


subroutine open 


ccccccccccccccoccccccccceccccccccccccccccccccccccccccccccccccccccccccccc 
cceececcceccceecccoocecccccccccccccccccccccccccccccccccccccccccccccccccc 
subroutine sopen(dummy) 
open(unlts1.filesOdataO/alpha") 
rewlndd) 

open(units2,files^dataO/Aindamental") 
reHlnd(2) 

open(unit33«flles'*data0/har2'*) 
rewind (3) 

open(units4«files'*data0/har3") 

rewlnd(4) 

open (units?»f iles"dataO/iaighar ") 
rewind(7) 

open(unitsS , files"dataO/hmag") 
rewind(8) 

open(units9,files"input") 

rewindO) 

return 

end 










ccoccccccccccccccccccececccccecocccecccccccccccccccccccccccccccccceccecc 
ccccccccccccccccccccccccccecccccecccecccceeeccccccccecccccceeccccccccccc 
c c 

c subroutine square c 

c c 

cccccccccccccecccccecceecccccecccccccceccccccccccccccccccccceecccccccccc 
cccccccccccccccccccccceccccceccccccccccccccccccccccccccccccccccccccccccc 
subroutine squareC f,fsquar,n ) 
dimension f(10) 
fsquarsO.O 
do 200 1=1,n 

fsquar3fsquar-t-f(l)*f(i) 

200 continue 
return 
end 
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oecccccccccccccccccccccccccccccccceccccccccccccccccccccccccccccccccccccc 
cccccccccccccccccccccccccceccccccccccccccccccccccccccccccccccccccccccccc 
c c 

c subroutine square c 

c c 

cccccccccccccecccccccccccccccccccccccecccccccccccccccccccccccccccccccccc 
occccccccccccccccccccccceccccccccecccccccccccccccccccccccccccccccccccccc 
subroutine square( f,fsquar,n ) 
dimension f(10) 
fsquarsO.O 
do 200 is1,n 

fsquarsfsquar+f(i)*f(i) 

200 continue 
return 
end 
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APPENDIX C 


FOURIER SPECTRUM PROGRAM 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c c 

c frequency spectrum c 

c c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c tor - Input data (this Is a real variable) 

c X - complex variable containing the Input data 

c before the fft subroutine call and 

c the complex frequency after the fft subroutine call 

c mag - magnitude of the complex frequency 

c n - number of data points (must be a power of 2) 

c Inv - used by fft 

c InvsO fft 

0 Invsl Inverse fft 

c maxi - used to decided highest frequency to be plotted 

c avemag - used to normalize data to 100} of the largest 

c harmonic 

dimension tor(1030) 
complex x(1030),cmplx 
real mag 
call sopen(O.O) 

pla3.1415926535 

rmssO.O 

read(9t*) n 

wrlte(6,*) ’ ns',n 

read(9,*) inv 

write(6,*) ' invs’.inv 

read(9f*) maxi 

wrlte(6,*) ' maxis',maxi 

read(9,*) iwrite 

write(6,*) ' iwritesiwrite 

do 100 ls1,n 

lf( inv .eq. 0 ) readd,2000) tor(i) 

lf( inv .eq. 0 ) x(l)scmplx( tor(i),0.0 ) 

lf( inv .eq. 1 ) read(4,2000) x(i) 

lf( Iwrlte .eq. 1 .and. Inv .eq. 0 ) wrlte(4,*) x(l) 

100 continue 

call fft( x,n,lnv ) 

do 200 Is1,maxi 
magscabs( x(l) ) 
pha8e3atan2( yy,xx ) 

avemagsamaxK cabs(x(1)),cabs(x(2)),cabs(x(3)),cabs(x(4)) ) 

rms3rms-HDag**2 

write(3»*) mag/avemag*100.0 













cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
cccccccccccccccccccccccccccccccccccccccccccccecccccccccccccccccccccccccc 
c c 

c subroutine close c 

c c 

cccccccecccccccccccccccccccccccccccccccccccccccccccccccecccccccccccccccc 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
subroutine sclose(duinmy) 
closed) 
close(2) 
close(3) 
close(4) 
close(7) 
close(8) 
close(9) 
return 
end 
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' 3 --Tl.' 


c 

o 

c 

c**The following fft subroutine Is taken almost verbatim 
c**from Ahmed K.* Rao K.R., Orthogonal Transforms for 
c**Dlgltal Signal Processing, Sprlnger^Verlag, 1975., pp79 
c 
c 
c 

subroutine fft(x,n,lnv) 

ceeeeeeeeeeeeeeeeeeeeeeeeeeefteeeeeeeeaeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeei 

c 

c this program Implements the fft algorithm to compute the discrete 
c fourler coefficients of a data sequence of n points 
c 

c calling sequence from the main program: 
c call fft(x,n,lnv) 
c n: number of data points 

c x: complex array containing the data sequence. In the end dft 

c coeffs. are returned In the array. Main program should 

c declare It as— complex x(1024) 

c Inv: flag for Inverse 

c InvsO for forward transform 

c Inval for Inverse transform 

c 

ceeeaeeeeeeaeeeaeeaeeieaaeftaeeaeeaaaaeaeaeeeaaaeeeeeeeeeeteeteeeeeieei 

complex x(1030),w,t,cmplx 
c 

c calculate the # of Iterations (log. n to the base 2) 

c 

ItersO 

Iremsn 

10 lremslrem/2 

If (Irem.eq.O) go to 20 
ltersiter+1 
go to 10 
20 continue 
slgns-1 

If (Inv.eq.l) signal 
nxp2an 

do 50 Ita1,Iter 
c 

o computation for each Iteration 

e nxp:number of points In a partition 

c nxp2:nxp/2 
nxpanxp2 
nxp2anxp/2 

wpwra3.141592/float(nxp2) 
do 40 ma1,nxp2 









calculate the multipier 

argsfloat(m-1)*wpwr 

wscmplx(cos(arg),sign*sln(arg)) 

do 40 mxpsnxp,n,nxp 

computation for each partition 

J1=mxp>nxp+m 

j2=jUnxp2 

t=x(j1)-x(J2) 

x(j1)=x(j1)+x(j2) 

40 x(J2)st»w 
50 continue 

unscramble the bit reversed dft coeffts. 

n2=n/2 

n1sn-1 

J=1 

do 65 isl.nl 
If(i.ge.J) go to 55 
t*x(j) 
x(J)sx(i) 

X(i)st 
55 ksn2 

60 if(k.ge.j) go to 65 
J»J-k 
k3k/2 

go to 60 
65 j=J+k 

if (lnv.eq.1) go to 75 
do 70 isl.n 
70 x(l)sx(i)/float(n) 

75 continue 
return 
end 












cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccoccccccccc 
c c 

c subroutine open c 

c c 

cccccccccccccccccccccccccccccccccccccccccccccecccccccccccccccccccccccccc 
cccccccccccccccccccccccccccccceccccccccccccccccccccccccccccccccccccccccc 
subroutine sopen(dummy) 
open<unlta1,flles"torque'*) 
rewlndd) 

open(unlts2,flles^dataO/phase") 
rewind(2) 

open(unlts3,fiie=*'dataO/Percent") 
rewind(3) 

open(unlts4,flles'»file1") 
rewind(4) 

open (unit a?, f 1 lea '•data0/mlac2'*) 
rewind(7) 

open(unlta8,fllea''data0/mlsc3”) 
rewind(8) 

open(unlta9,filea"Input") 

rewind(9) 

return 

end 








APPENDIX D 


PERFORMANCE PROGRAM 


cccccccccccccccccccccccccccccccccccccccecccccccccccccccocccccccccccceccc 

cccccccccccccccccccccccccccccccccccccccccccccecocccccccccccccccccceccccc 


c c 
c program to calculate developed power- c 
c c 
c speed points for brushless dc motor c 
c c 


cceccceccccccccccccccccccceccccccccccccccccccccceccccccccccccccceccccccc 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

c 

c rpm - speed of dc motor 

c wm - angular speed of dc motor 

c p - number of poles 

c t - time 

c kO - number of Intervals between print outs 

c h - interval of integration 

c w - electrical angular frequency of pan 

c ws - electrical angular frequency of pmg 

c alpha - commutation advance angle (degrees) 

c ra - stator resistance 

c xla - stator leakage inductance 

c rO - choke coil resistance 

e xlo - choke coll inductance 



xk - motor constant 

xkg - generator constant 

vd - Inverter input voltage 

X - array of stator phase currents 

tp - period 

ang - electrical angular displacement 

xlamda - electrical angular displacement (pmg) 

m - number of phases 

pd - total developed power 

td - total developed torque 

pdl - phase 1 developed power 

pd2 - phase 2 developed power 

pd3 - phase 3 developed power 

e(1) - excitation voltage of phase 1 

e(2) - excitation voltage of phase 2 

e(3) • excitation voltage of phase 3 


dimension f(6),v(6)»x(6),xx(6),a(6),r(6),e(6),alphs(2) 

common rr»rd,rO,ra,r(6),piyXla»xlO»e(6),vm,vdc,ivdon,imodon 

real iav.irms 

dummysO.O 

call popen(dummy) 

pi*3.1415926536 

eonvapi/180.0 

read(9»*) istop 
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write(6,*) ' istops'•istop 

read(9,*) alpha 

write(6,*) ' alphas',alpha 

read(9,*) ( a(i).lsl,5 ) 

write(6,«) ' as',( a(i),isi,5 ) 

read(9,*) h 

wrlte(6,*) ' hs'.h 

read(9,*) rO 

write(6,») » rOs'.rO 

read(9,») xlO 

write(6,•) » xlOs'.xlO 

read<9,*) ra 

write(6,'») ' ras'.ra 

read(9.*) xla 

write(6,») * xlas'.xla 

read(9,*) rd 

write(6,») » rds'.rd 

read(9«*) rpmO 

write(6,*) * rpinOs',rpmO 

read(9,») xk 

write(6,») » xks«,xk 

read(9,») xkg 

write(6,*) ' xkgs'.xkg 

read(9»*> pmax 

write<6,») ' pfflaxs'.pnax 

read(9f*) pmin 

write<6,») » pmins'.pmin 

read(9f*) vdc 

write(6,*) ' vdes',vdc 

read(9f*) points 

write(6,*) ' points**,points 

read(9,*) ivdon 

writeC6,») » ivdonsivdon 

read(9»*) ioodon 

write(6,*) ' imodons',iinodon 

read(9f*) ichon 

write(6,») ' ichons',ichon 

read(9**) iwgon 

write(6,») ' iwgonsiwgon 

read(9,») ( x(i),isi ,3 ) 

write(6,») » xs*,( x(i),isi,3 ) 

read(9,*) iarea 

write(6,*) ' iareas',iarea 

read(9«*) xwg 

write(6,*) ' xwgs<,xwg 

read(9»*) ifull 

write(6,») • ifulls',ifull 

do 4 Js1,5 

lf( ioodon .eq. 1 ) a(J)s a(J)»pi/l80.0 
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'i L'- *r *A't".<T^-^.'':^'.'r.^'.v.'7.'^.^.‘'T.»r:T’.*r^ 


1 


0 ) a(j)=pi/2.0 
0 ) a(1)s30.00»pi/180.0 


) wg=7950»pi 
) v(gsxwg*v)m 


lf( Imodon .eq. 
if( imodon .eq. 
continue 
ms3 

pfw=(rpmO/45000.0)••2«330.0 
wmsrpm0*pi/30.0 
lf( Iwgon .eq. 1 
if( iwgon .eq. 0 
wmsabs(wm) 
ps4.0 
t=0.0 

wsp/2.0*wm 

wss2.0*wg 

vm=sqrt(3.0)*xkg*wg 
wrlte(6,») ’ vms'.vm 
alphasalpha*conv 

areas0.5-i>0.5*( pi/3.0-a(1 )<t>a(2)-a(3)->-a(4)-a(5) )/( pi/3.0 ) 


0 .and. iarea .eq. 1 ) vdcsvde/area 
eq^ 1 ) alph3(1)saco8( (vde«sqrt(2.0))/(1.0»vm»1.35) ) 
eq. 1 ) alphs(2)saeoa( (vdc*8qrt(2.0))/(2.0*vm*1.35) ) 
eq. 0 ) alphs(1)sacos( (vdc»aqrt(2.0))/(0.5»vm»1.35) ) 
eq. 0 ) alphs(2)8aeos( (vdc«sqrt(2.0))/(1.0»vm*1.35) ) 
' alphs(1)s',alphs(1)/conv 
' alphs(2)s',alphs(2)/eonv 
' areas*.area 
' f*»,w/(2.0»pi) 

* fsa',ws/(2.0»pi) 

' mag of e1s',xk*w/2.0 


lf( imodon 
if( ifull 
lf( iftill 
if( ifull 
if( ifull 
wrlte(6,*) 
write(6,») 
write(6,») 
write(6,») 
write(6,*) 
wrlte(6,») 
xxCDsO.O 
xx(2)s0.0 
xx(3)-0.0 
pdsO.O 
pinsO.O 
iavsO.O 
irmasO.O 
elavesO.O 
vdavesO.O 
eIrmssO.O 
sumsO.O 
tps2.0*pi/w 
tminspmln*tp 
tmaxspmax*tp 
tmxmnstmax-tmin 
soales30.0/pointa 

kOsintC tmxmn/tp*8cale/ab8(rpmO)/h ) 

kountsO 

index*1 

rewind(9) 

read(9**) ietop 

if( i8top .eq. 1 ) go to 305 







kountskount-fl 

angsv*t 

Iprin1 Sint(ang/2.0/pi) 

angsang-float(iprin1)*2.0*pi 

a(1)sxk*w*2.0/p*3in(ang-alpha) 

e(2)sxk»w»2.0/p*sin(ang-alpha-2.0/3.0*pi) 

e (3 )sxk*w*2.0/p*ain (ang-alpha-i>2.0/3.0*pi) 

wst3wa*t0 

if( wst .ge. 2.0*pi ) tOsO.O 
wstsws*tO 

if( ifull .eq. 0 ) call half( MSt.xlamda ) 
if( ifull .eq. 1 ) call fulK wst.xlamda ) 
call fv( a,ang«v,vd,alphs,xlaiDda ) 
call ff( indextt,f,x,xl,ang,del»v ) 
if( ichon .eq. 1 ) call check( m,x,f ) 
call derk( m,Xtf,t,tO,h(lndex ) 
if( index .eq. 2 ) go to 200 
if( t .It. tmin ) go to 100 
f(1)s( x<1)-xx(l) )/h 
f(2)s( x(2)-xx(2) )/h 
f(3)s< x(3)-xx(3) )/h 
do 86 kks1,3 
xx(kk)sx(kk) 
continue 

call ffl( fi,f,ang ) 

call fvap( wst,vap,alphas ) 

v1psvd-fi»xl0/2,0-xi«r0-x(1)»(r(1)-ra) 

vscrsvap-vlp 

if( ang .It. 3.0/3.0*pi .and. wst .It. 2.0/3.0«pi ) 
vcsrsx(l)»(r(1)-ra) 

if( ang .ge. 3.0/3.0*pi .and. wst .ge. 3.0/3.0*pi .and. 

wst .It. 4,0/3.0»pl ) vscrsx(l)*(r(1)-ra) 
lf( t .ge.tmax ) go to 300 
if( kount .It. kO ) go to 100 
sumssum^'l. 0 
pd1sx(1)*e(1) 
pd2sx(2)»e(2) 

pd3sx(3)»e(3) 

tds( pd1-t>pd2<*-pd3 )/wn 

pdspd-*-( pd1-*>pd2-*-pd3 ) 

pinspin4>( v(1)»x<1)+v<2)»x(2)-*-v(3)*x(3) ) 

lrmssirms-t.x(1)*»2 

iavsiav4-abs( x(1) ) 

e 1 nasse 1 rms4>e (1) **2 

e1avese1ave<»>abs( e(1) ) 

vdavesvdave-t>vd 

h8in6shsln6-f td*sin (6.0*ang) 

hoo86shco86-«>td*oos (6.0*ang) 

h8in12sh8in12<*>td*8in(12.0*ang) 

hoo812shoo812-ftd*oo8 (12.0*ang) 







hsin18shsin18-t-td*sln(18.0«ang) 

hcos 18shcos 18-i>td*cos (18.0*ang) 

hsln24shsin24+td«sln(24.0«ang) 

hcos24shcos244-td*cos (24.0*ang) 

txst-toxmn 

txxst»w»180,0/pi 

write(1,2000) x(1) 

wrlte(2,3000) txx 

write(3,2000) td 

write(4,2000) v(1) 

write(7,2000) vd 

wrlte(8,2000) x(1)+x(2)+x(3) 

kountsO 

lf( t .It. tmax ) go to 100 
305 pdspd/sum 
pinspin/sum 
irmaslrma/sum 
IrmsssqrtC irna ) 
iavsiav/aum 
e1rmsxe1rms/sum 
e1rms3sqrt( eirma ) 
elavexelave/suffl 
vdavesvdave/suffl 
haln6«haln6/sum 
hcos6shco86/suiB 
h6shsin6”2+hcos6”2 
h6saqrt( h6 )/6.0 
hsln12shain12/su« 
hcos12shcos12/sum 
h12*hsin12««2+hcos12»«2 
h12asqrt( h12 )/12.0 
hslnl8shsin18/sum 
hoos18shcos18/sum 
hI83hslnl8**2«hcos18**2 
h18ssqrt( h18 )/18.0 
hsln24shsin24/sum 
heos24shcos24/sum 
h24shsin24**2i4icos24**2 
h24asqrt( h24 )/24.0 
600 tsav3(pd-pfw)/wm 

write(6,*) ’ tsavs*,tsav 
effx(pd-pfw)/pin*100.0 
if< eff .It. 100.0 ) go to 710 
ef f spin/(pd4>pfw) * 100.0 
710 wrlte(6,*) ' pins’,pin 

wrlte(6,*) ' effs’,eff 
write(6,*) ' lavs',lav 
wrlte(6,*) ' Irmas',irms 
wrlte(6,*) ' elaves',slave 
write(6,*) ' vdaves',vdave 
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Hrite(6,») • h6=',h6 
write(6,») • h12s*,h12 
write(6,») • h18s',h18 
write(6,») • h24s»,h24 
1000 format( 5(e12.4,5x) ) 

2000 formate e10.4 ) 

3000 formate f8.3 ) 

call pcloseedummy) 

stop 

end 
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cccccecocccoceccccccceeccccccccccccccecccccccccccccccccccccccccccccccccc 
ccccccccccceccccccccccccccccccccccccccccccccccccccccccccccccccccccccccec 
c c 

c subroutine check c 

c c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
subroutine check( m,x,f ) 
dimension x(6),f(6) 

common rr,rd,rOyra.r(6),pl,xla,xlO,e(6),vm,vdc,lvdon,lmodon 
do 10 lsl,m-1 

slgnfsabs( f(l) )/f(l) 
slgnxsabsC x(l) )/x(l) 

lf( abs( f(l) ) .gt. 1.0e08 ) f(l)sslgnf*1.0e08 
lf( abs( x(l) ) .gt. 1.0e06 ) x(l)sslgnx*1.0e06 
10 continue 

f(3)«-( f(1)+f(2) ) 

return 

end 






cceccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

cccccccccccccoccccoccccccccccccccccccecccccccccccccccccccccccccccccccccc 


i 


c subroutine pclose c 

c c 

cccecccccccccccccccccccccccccccccccccccccccccccccccccccccccccccocccccccc 
cccccccccccccccccceecccccccccccccccccccceeeccccccccccccccccccccccccccccc 
subroutine pclose(dummy) 
close(1) 
close(2) 
close(3) 
close(4) 
close(7) 
close(8) 
close(9) 
return 
end 









cccccccccccccccccccccccccccccccccccccecccccccccccccccccccccccccccccccccc 

cccccccccccccccccccccccccceccccccccccccceccccccccccccccccccccccccccccccc 


c c 

c subroutine derk c 

c c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c derk Is a fourth-order, fixed increment runge-kutta 

c integration routine 

c 

c m 

c X 

c f 

c h 

c t 

c index 

c 1. if 

c 2. if 

c 

subroutine derk( m,x,f,t,tO,h,index ) 
dimension x(6),f(6),q(400) 
lf( index .eq. 2 .) go to 19 

18 kxxsO 
lndexs2 
do 35 is1,m 
jsi+300 

35 q(J)sx(l) 

19 kxxskxx-i-1 
go to (1,2,3t4),kxx 

1 do 5 is1,m 
q(i)sf(i)»h 

5 x(i)sx(i)+q(i)/2.00 
tst+h/2.00 
tOstO-fh/2.00 
return 

2 do 6 is1,m 
Jsi4.100 
ks 14-300 
q(J)sf(i)»h 

6 x(l)sq(k)+q(J)/2.00 
return 

3 do 7 i«1,m 
Jsi-^00 
ksi-i-300 
q(J)*f(i)»h 

7 x(i)*q(k)+q(J) 
t>t4>h/2.00 
t0st044)/2.00 
return 


- number of simultaneous differential equations 

- array of dependent variables 

- array of derivatives of independent varibles 

- Increment 

- time 

- Indicator 

exit routine with index=1, solution found at=t+h 
exit routine with lndexs2, go back to re-evaulate 





do 8 is1,m 
jsi+100 
ksi-i-200 
l=i+300 

x(i)sq(l)+( ( q(l)45.00»q(k)+2.00»q(k)+f(i)»h )/6.00 ) 

indexsl 

return 

end 








cccccccccccccccccccccccccccocccccccecccccccccccccccccccccccccccccccccccc 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


cccccccccccccccccccccccccccceccecceeeccccccccccccccccccccecccceccccccccc 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccecc 

subroutine ff( index,t,f,x,xl,ang,del,v ) 
dimension f(6),v(6),x(6) 

common rr,rd,r0,ra,r(6),pi,xla,xl0,e(6),vm,vdc,ivdon,ifflodon 
rrsrd+ra 

if( index .eq. 1 ) r(1)srr 
if( index .eq. 1 ) r(2Urr 
if( index .eq. 1 ) r(3)=rr 

if( ang .ge. 0.0 .and. ang .It. 1.0/3.0«pi ) goto 10 

if( ang .ge. 1.0/3.0*pi .and. ang .It. 2.0/3.0«pi ) goto 20 

if( ang .ge. 2.0/3.0«pi .and. ang .It. 3.0/3.0»pi ) goto 30 

if( ang .ge. 3.0/3.0«pi .and. ang .It. 4.0/3.0«pi ) goto 40 

if( ang .ge. 4.0/3.0«pi .and. ang .It. 5.0/3.0*pi ) goto 50 

if( ang .ge. 5.0/3.0»pi .and. ang .It. 6.0/3.0*pi ) goto 60 

10 nisi 

n283 
n3s2 

if( icyele .eq. 1 ) iflagsO 
itestalO 
icyclesO 
go to 100 
20 nixl 

n2a2 
n3a3 

itestx20 
go to 200 
30 n1x2 

n2xi 

n3*3 
itestx30 
go to 100 
(0 n1x2 

n2x3 

n3«1 

itestx40 
go to 200 
>0 n1x3 

n2x2 
n3«1 

itestxSO 
go to 100 









n2sl 
n3s2 
itests60 
Icyclesl 
go to 200 

100 lf( x(n2) .It. ^.001 .and. index .eq. 1 ) iflagsitest 

lf( x(n1) .It. -0.001 .and. index .eq. 1 ) r(n1 )srr-t-1000.0*abs( x(n1) ) 

if( x(n3) .gt. 0.001 .and. index .eq. 1 ) r(n3)srr+1000.0«abs( x(n3) ) 

if( iflag .eq. itest .and. index .eq. 1 ) r(n2)srr-t-1000.0*abs( x(n2) ) 

do 150 is1,3 

if( r(i) .gt. 1000.0 ) r(i)=1000.0 
150 continue 

xlsxla 

f(n1)s( -x(n1)*r(n1)4-v(n1)-e(n1) )/xl 
f(n2)s( -x(n2)«r(n2)+v(n2)-e(n2) )/xl 
f(n3)s( -x(n3)»r(n3)+v(n3)-e(n3) )/xl 
return 

200 if( x(n2) .gt. 0.001 .and. index .eq. 1 ) iflagsitest 

if( x(n3) .gt. 0.001 .and. index .eq. 1 ) r(n3)»rr<«>1000.0*abs( x(n3> ) 

if( x(n1) .It. -0.001 .and. index .eq. 1 ) r{n1)xrr+1000.0«abs( x(n1> ) 

if( iflag .eq. itest .and. index .eq. 1 ) r(n2)srr4>1000.0*abs( x(n2) } 

do 250 ix1,3 

if( r(i) .gt. 1000.0 ) r(i)*1000.0 
250 continue 

xlsxla 

f(n1)=( -x(n1)»r(n1)+v(n1)-e(n1) )/xl 
f(n2)x( -x(n2)*r(n2)4-v(n2)-e(n2) }/xl 
f(n3)*( -x(n3)*r(n3)+v(n3)-e(n3) )/xl 

return 

end 







ccccocccccoocccccccccececccccccccccocccccccccccccccccccccccccccceccccccc 
cccceeccccccccceececcccccccceccecccccccccccccccccccccccccccccccecccccccc 
c c 

c subroutine fl c 

c c 

cccccccccccccccccccccccccccceccccccccccccccceccccccccccccccccccccccccccc 
ccccccccccoeececcccccccccccccecccccccccccccccccccccccccccccccccccccccccc 

subroutine ffl( fl,f,ang ) 
dimension f(6) 

common rr,rd,r0,ra,r(6),pl,xlatxl0,e(6),vm,vdc,lvdon,lmodon 
lf( ang .ge. 0.0/3.0»pl .and. ang .It. 1.0/3.0«pl ) fl=-f(2) 

lf( ang .ge. 1.0/3.0«pl .and. ang .It. 2.0/3.0«pl ) fls f(1) 

lf( ang .ge. 2.0/3.0*pl .and. ang .It. 3.0/3.0*pl ) fls-f(3) 

lf( ang .ge. 3.0/3.0«pl .and. ang .It. 4.0/3.0«pl ) fls f(2) 

lf( ang .ge. 4.0/3.0»pl .and. ang .It. 5.0/3.0«pl ) fl=-f(1) 

lf( ang .ge. 5.0/3.0«pl .and. ang .It. 6.0/3.0«pl ) fls f(3) 

return 
end 










ccccccccccecccccoocccccccccccccccccccecccccccccccccccccccccccccccccccccc 
cccccccccccccccccccccccecccccccccccccccecccccccccccccccccccccccccccccccc 
c c 

c subroutine full-bridge c 

c c 

ccccceccccccccccccececeecccccoceceecoccccccccccccccccccccccccccccccccccc 
cecccccccccccccccccccccccccccccccccccccccccccccccccceeeocceccecccceccccc 

subroutine fulK ust.xlamda ) 

common rr,rd,rO«ra,r(6),pi,xla,xlO,e(6),vm,vdc,lvdon,lmodon 


lf( wst .ge. 0.0 
wst-0.0 

if( wst .ge. 1.0«pi/3.0 
wst-1.0«pi/3.0 
if( wst .ge. 2.0*pi/3.0 
wst-2.0«pl/3.0 
if( wst .ge. 3.0«pl/3.0 
wst-3.0»pi/3.0 
lf( wst .ge. 4.0*pl/3.0 
wst-4.0»pi/3.0 
lf( wst .ge. 5.0«pl/3.0 
wst-5.0»pi/3.0 
return 
end 


.and. wst .It. 1.0*pl/3.0 ) xlamdas 
and. wst .It. 2.0*pl/3.0 ) xlamdas 
and. wst .It. 3.0*pl/3.0 ) xlamdas 
and. wst .It. 4.0*pi/3.0 ) xlamdas 
and. wst .It. 5.0*pi/3.0 ) xlamdas 
and. wst .It. 6.0*pi/3.0 ) xlamdas 







cccccccccccccccccccccccccccccccceccccccccccccccccccccccceccccccccccccccc 
cceccccccccccccccccccccccccccccccccccccccccccceccecccccccccccccccccccccc 
c c 
c subroutine half-bridge c 
c c 
cccccccccccccceeeecccecccccccccccccccccccccccceccccceececccccccccccccccc 
cccccccececccccccecccccccccccoeeeeececcccccccccccccccccccccccccccccccccc 


subroutine half( wst.xlamda ) 

conunon rr,rd,r0,ra,r(6),pl»xla,xl0,e(6),vm,vdc,lvdon,iraodon 
if( wst .ge. 0.0 .and. wst .It. 2.0*pi/3.0 ) xlamdas 

c wst-0.0 

if( wst .ge. 2.0*pi/3.0 .and. wst .It. 4.0*pi/3.0 ) xlamda= 
0 wst-2.0*pl/3.0 

lf( wst .ge. 4.0*pl/3.0 .and. wst .It. 6.0*pi/3.0 ) xlainda= 
c wst-4.0*pi/3.0 

return 
end 
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open(units3tfiles"dataO/Torque") 
rewind(3) 

open(unlt84,file«"dataO/Van'*) 
rewind(4) 

open(units?,files-dataO/Vd") 
rewind(7) 

open(units8,flles«dataO/In") 
rewind(8) 

open(units9,files"input*) 

rewind(9) 

return 

end 
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cceeccccccccccecceeccoecccccecececccccceeececccccccccccccccccccccccccccc 
ccccoccocceccccccccccccccccccccccccccccccccccccccceccccccccccccccccccccc 
c c 

c subroutine vd c 

c c 

cceeccceocccoccccececoceecccceeoeeeeeeoecccccccccecccecccccecccccccccccc 
cccccccccccceccccccccccccococccccccccccccccccccccccccccccccccccccccccccc 
subroutine fv( a,ang,v,vd,alphs«xlafflda ) 
dimension a(6),vU),alphs(2) 

common rr,rd,r0,ra,r(6),pi,xla«xl0ye(6),vm,vdc»ivdon,imodon 
betasxlamda-t-alphs( 1) 

if( ivdon .eq. 1 .and. ifull .eq. 1 ) vd=vm*sin( beta<fpi/3.0 ) 

if( ivdon .eq. 1 .and. ifull .eq. 0 ) vdsvm*sin( beta-fpi/6.0 ) 

if( ivdon .eq. 0 ) vdsvde 

ahift=30.0/180.0»pi 

thetasshlft-fang 

call svxnC vx,theta,a ) 

vansvx 

thetasshif t4>ang-2.0/3.0*pi 
call svxn( vx,theta.a ) 
vbnavx 

theta«shift4>ang-«>2.0/3.0*pi 
call svxn( vx,theta,a ) 
vonsvx 

v(1)»vd*van 

v(2)avd*vbn 

v(3)»vd*vcn 

return 

end 
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subroutine fvap( W8t,vap,alps ) 

eonmon rr,rd,rO,ra,r(6),pl»xla,xlO,a(6),vm,vdc,ivdon,iiBodon 
vaaaO.O 


vabs vm*sln( wst4>1.0/3.0*pi-»-alp8 ) 
vaes-vni*sln( wst4>5.0/3.0*pl<falp8 ) 


if( 

ust 

.ge. 

0.0/3.0«pi 

.and. 

wst 

.It. 

1.0/3.0»pi 

) vap« 

vab 

if( 

wst 

.ga. 

1.0/3.0«pi 

.and. 

wst 

.It. 

2.0/3.0«pl 

) vapx 

vac 

if( 

wst 

.ge. 

2.0/3.0*pi 

.and. 

wst 

.It. 

3.0/3.0»pi 

) vaps 

vac 

if( 

wst 

.ga. 

3.0/3.0»pl 

.and. 

wst 

.It. 

4.0/3.0»pl 

) vapx 

vaa 

if( 

wst 

.ga. 

4.0/3.0«pl 

.and. 

wst 

.It. 

5.0/3.0«pi 

) vapx 

vaa 

if( 

wst 

.ga. 

5.0/3.0*pl 

.and. 

wst 

.It. 

6.0/3.0»pl 

) vapx 

vab 


raturn 

and 
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cccccccccccccccceceecccccccccceccccccccccccccccccecccccccccccccccccccccc 

ccccccecccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccec 


c subroutine vxn 

c 


c 

c 

c 


ccccccccoccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
subroutine svxnC vx,th,a ) 
dimension a(6) 

common rr,rd,rO,ra,r(6),pi,xla,xlO,e(6),vm,vdc,ivdon,imodon 

vxsO.O 

xasO.O 

xbspi 

xcs2.0*pl 

a(6)apl/2.0 

thsth+2.0»pi 

th=th-float( int(th/xc) )*xc 

lf( th .ge. xa-t-ad) .and. th .It. xa-t>a(2) ) vxsi.O 

if( th .ge. xa4-a(3) .and. th .It. xa-fa(4) ) vxsi.O 

lf( th .ge. xa-«>a(5) .and. th .It. xa<t>a(6) ) vxsi.O 

if( th .ge. xb-a(6) .and. th .It. xb-a(5) ) vxsi.O 

if( th .ge. xb>a(4) .and. th .It. xb-a(3) ) vxsi.O 

if( th .ge. xb-a(2) .and. th .It. xb-a(1) ) vxsi.O 

lf( th .ge. xb4>a(1) .and. th .It. xb-«-a(2) ) vxs.1.0 

lf( th .ge. xb<^a(3) .and. th .It. xb-*>a(4) ) vxs-l.O 

if( th .ge. xb4>a(5) .and. th .It. xb4-a(6) } vxs>1.0 

if( th .ge. xc<-a(6) .and. th .It. xe-a(5) ) vxs-1.0 

if( th .ge. xc-a(4) .and. th .It. xe>a(3) ) vxs-1.0 

if( th .ge. xc-a(2) .and. th .It. xc>a(1) ) vxs.1.0 


return 

end 







